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ABSTRACT 

A method for improving the accuracy of hydrodynamical codes that use a moving Voronoi 
mesh is described. Our scheme is based on a new regularization scheme that constrains the 
mesh to be centroidal to high precision while still allowing the cells to move approximately 
with the local fluid velocity, thereby retaining the quasi-Lagrangian nature of the approach. 
Our regularization technique significantly reduces mesh noise that is attributed to changes in 
mesh topology and deviations from mesh regularity. We demonstrate the advantages of our 
method on various test problems, and note in particular improvements obtained in handling 
shear instabilities, mixing, and in angular momentum conservation. Calculations of adiabatic 
jets in which shear excites Kelvin Helmholtz instability show reduction of mesh noise and 
entropy generation. In contrast, simulations of the collapse and formation of an isolated disc 
galaxy are nearly unaffected, showing that numerical errors due to the choice of regularization 
do not impact the outcome in this case. 
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1 INTRODUCTION 

New computational methods for fluid dynamics that employ a mov¬ 
ing Voronoi mesh approach have been developed in the past several 
years to simulate astrophysical and cosmological systems. Codes 
based on this method include AREPO (Springel 2010), which 
has recently been used to run the 12 billion resolution element 
state-of-the-art “Illustris” cosmological simulation (Vogelsberger 
et al. 2014a,b; Genel et al. 2014), TESS (Duffell & MacEadyen 
201 1), adapted for relativistic hydrodynamics, the moving mesh al¬ 
gorithm presented in Gaburov, Johansen & Levin (2012) for the 
simulation of magnetically levitating accretion disks around super- 
massive black holes, and the RICH code (Steinberg et al. 2015). 
The moving Voronoi framework has also been extended to finite- 
element techniques (Mocz et al. 2014) and constrained transport 
approaches for magnetohydrodynamics (MHD) to strictly maintain 
the divergence-free nature of the magnetic held (Mocz, Vogels¬ 
berger & Hernquist 2014). Related mesh-free methods have also 
been developed by Hopkins (2014), which use Riemann solvers 
acting over volume “overlaps”, and have similar advantages to the 
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moving mesh approach, and have particularly good angular mo¬ 
mentum conservation but at the cost of enhanced noise. 

Moving mesh codes have found success due to the numerical 
advantages they provide because of their quasi-Lagrangian nature, 
Galilean-invariant formulation (for non-relativistic fluids), limited 
advection errors, preservation of contact discontinuities, contin¬ 
uous spatial adaptability, and ability to accurately resolve insta¬ 
bilities. Notably, the moving mesh method provides advantages 
over previous smoothed particle hydrodynamics (SPH) approaches, 
which can suppress entropy generation by mixing, underestimate 
vorticity generation in curved shocks, prevent efficient gas stripping 
from infalling substructures (Sijacki et al. 2012), and have rela¬ 
tively poor convergence properties (e.g. Zhu, Hernquist & Li 2015); 
and also over adaptive mesh refinement (AMR) codes, which can 
have large advection errors due to numerical diffusion in the pres¬ 
ence of large supersonic bulk velocities. 

However, moving Voronoi mesh codes can be affected by 
noise on small spatial scales due to the mesh motion (Bauer & 
Springel 2012; Hopkins 2014). As the mesh evolves, the orien¬ 
tation and sizes of the faces and the topology of the connections 
between cells change, introducing small truncation errors in the so¬ 
lution, which can lead to artificial additional power on the smallest 
scales (Bauer & Springel 2012) and numerical secondary instabil- 
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ities in shear flows. The errors are largest in the presence of shear, 
when the mesh cells move past each other and change face bound¬ 
aries rapidly, because the standard moving mesh formulations as¬ 
sume that the face orientations do not evolve signiflcantly in one 
timestep or may not adequately capture rapid changes. If the mesh 
generating points strictly move with the fluid velocity, the mesh 
noise can be quite large because the Voronoi cells can deform into 
highly irregular shapes. 

Mesh regularization, i.e., applying a small correction to the 
mesh generating points of the cells in a way that keeps the cells 
“well-behaved” (round and centroidal) provides a handle on this 
error. The mesh generating points travel with the fluid velocity plus 
some small correction to keep the mesh approximately regular and 
centroidal. The approach taken by Springel (2010) adds a small cor¬ 
rection in the direction of the center-of-mass at the beginning of the 
time step (or sub time step in the case of higher-order time integra¬ 
tion). Vogelsberger et al. (2012) modify this approach for cosmo¬ 
logical simulations by considering the maximum opening angle of a 
face at the mesh generating point and the density gradient direction 
to set the correction. The correction term scales as the sound-speed 
inside the cell, which ensures a timestep independent formulation. 
We will refer to this class of regularization methods as sound-speed 
regularization (SSR). 

Despite keeping the cells approximately regular, these tech¬ 
niques do not eliminate mesh noise entirely because there are sev¬ 
eral small sources of error present in the formulation. First, the mo¬ 
tion of the mesh generating points (and hence the change in the 
orientation of the faces) is not smooth. The direction of the correc¬ 
tion only takes into account the geometry at the beginning of the 
time step. After the mesh advances to the next timestep, a small 
offset between the mesh generating point and the center-of-mass 
remains, which may point in an uncorrelated direction from that at 
the beginning of the time step. The cell then obtains an small kick 
in that uncorrelated direction, adding a small level of random noise 
to the smooth mesh deformation. Second, a small offset exists in the 
center-of-mass and mesh generating point. This causes the cells to 
have some spin about their center of mass, and also introduces er¬ 
rors in the second-order estimates of cell gradients (Springel 2010), 
which assumes the two points coincide. This later issue has recently 
been resolved in Pakmor et al. (2015), where an improved least- 
squares-flt gradient estimator has been developed that can achieve 
second-order accuracy for smooth flows. 

Recently, Duffell & MacFadyen (2014) proposed a smoothing 
of the velocities of the mesh generating points to alleviate the prob¬ 
lem of mesh noise. Their scheme exhibits a reduction in: grid noise, 
numerical secondary instabilities, and artiflcial power in the power- 
spectrum of the solution held on the smallest spatial scales (al¬ 
though not at the level of flxed grid codes). The smoothing method 
overcomes the issue of the uncorrelated randomly directed correc¬ 
tions that violate a fully smooth mesh motion. However, setting 
the cell vertex velocities to the fluid velocity and smoothing by the 
neighbours does not enforce a strongly centroidal mesh (only ap¬ 
proximately centroidal; as do the other regularization methods in 
the literature); deviations from a centroidal mesh introduce errors 
in e.g. Green-Gauss gradient estimates. Moreover, in the presence 
of shear flows, the cell velocities in this approach will be smoothed 
out, reducing the code to a static mesh code, and thus losing the 
desired properties of a moving mesh code such as advecting the 
solution at high precision. 

We have developed a new regularization scheme to address 
the smoothness and strong-centroidality issues, which is presented 
here. We refer to the method as strongly-centroidal Lloyd regu¬ 


larization (LR). Our approach evolves the mesh in a smooth, cen¬ 
troidal manner. Unlike the previous techniques, it accounts for 
where the center-of-mass of a cell travels to at the end of a time 
step in order to ensure the centroidal property, making it in essence 
a forward-predicting, iterative Lloyd’s algorithm. The concept is 
straightforward. The velocities of mesh generating points are ini¬ 
tialized to the values of the local fluid velocity, and a few iterations 
of corrections are applied to modify this velocity so that the cell 
moves to the location of its center-of-mass at the end of the time 
step. This results in a smoothly evolving centroidal Voronoi mesh 
that moves approximately with the fluid velocity. The scheme may 
optionally be modifled to a weighted centroidal scheme (weighted 
by the fluid density) for simulations with collapsing structure to 
gain automatic increased resolution in regions of high density. We 
demonstrate signiflcant reduction of mesh noise with the new reg¬ 
ularization scheme. 

Developing reliable numerical codes that offer general adapt¬ 
ability to a huge spatial and temporal dynamic range is an ongoing 
challenge. Recently Pakmor et al. (2015) have improved the accu¬ 
racy and angular momentum conservation for moving-mesh codes 
by developing new time integration scheme and spatial gradient es¬ 
timates in the AREPO code, signiflcantly improving the precision 
of the code in smooth test problems. In that work, the MUSCL- 
Hancock (MH) scheme (Toro 1999) is replaced by a second-order, 
time symmetric Runge-Kutta (RK) integrator via Heun’s method. 
Additionally, the Voronoi-improved Green-Gauss (GG) gradient 
estimates are replaced by a least-squares gradient estimator (LSF). 
Our regularization complements these recent improvements for 
smooth flows by improving the accuracy of moving mesh algo¬ 
rithms when shear is present. 

We describe the method and its computationally efficient im¬ 
plementation in Section 2. Several numerical tests are performed 
to show the advantages of our approach, the results of which are 
presented in Section 3. We offer concluding remarks in Section 4. 


2 STRONGLY-CENTROIAL LLOYD REGULARIZATION 

The regularization scheme involves initially setting the velocity of 
the mesh generating point of each cell i equal to the local fluid 
velocity : 

Wi,o = Vi (1) 

Then, in each successive iteration J, the predicted locations of the 
centers-of-mass of the cells Ci,j are computed at the end of the 
timestep At, and the mesh generating point velocities are updated 
to: 

Wi,j = (Ci,j - Fi) /At (2) 

where is the location of cell i at the beginning of the timestep. 

This iteration is essentially a Lloyd’s algorithm for converging 
an arbitrary Voronoi diagram to a centroidal one. We take 5 itera¬ 
tions in the simulations shown, as the method quickly produces a 
strongly centroidal mesh. We have explored this free parameter and 
found that 2 iterations is sufficient to produce quantitatively very 
similar results. 

Optionally, a limiter on the magnitude of the correction to the 
velocity may be added to explicitly enforce the motion of the mesh 
generating particles to be close to Lagrangian. We propose the fol¬ 
lowing, Galilean-invariant limiter, based on the local sound speed 
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Figure 1. Comparison of KHI with the new (LR) and previous (SSR) regularization schemes at low and high resolution at f = 2.0. LR eliminates secondary 
numerical instabilities, which are present in the simulations that use SSR, even with the recent RK and LSF improvements. For comparison, runs on a static 
grid are also shown, which also do not have secondary numerical instabilities. 


Cs,i of the cell. 


Ci^j ^ (r^ + Atwi,o) + (ci,j - (r^ + Atwi,o)) 


X min 



fcs^jAt _ 

(r^ + AtWi,o)|| 


( 3 ) 


That is, we restrict the motion of the cell to remain within a circle 
of radius /cs,i At, (0 < / < 1) around the first predicted location 
of the mesh vertex (which is just predicted by the fiuid velocity). 

The scheme may also be modified so that the centers-of-mass 
Ci^j are calculated by weighting with the density field, to bias the 
movement of the cells to regions of high density. Alternatively, for 
a simpler implementation with the same effect, the initial velocities 
in Equation 1 may be biased by a density gradient term as in Vo- 
gelsberger et al. (2012), to move cells together towards regions of 
collapse. 

An efficient implementation of the method does not require 
new mesh reconstruction in each iteration. The center-of-mass of 
cell i in the predictive step may be calculated by just using the 
points that are neighbours of the cell at the beginning of the 
timestep (kicking them forward in time to build the predictive 
Voronoi cell). This strategy avoids additional in-circle tests. In 
cases where the mesh topology changes, this is not strictly an exact 
calculation of the center-of-mass, but the error is negligible since a 
Voronoi diagram changes continuously. 


3 NUMERICAL TESTS 

3.1 Reduction of artificial secondary instabilities in the 
Kelvin Helmholtz instability 

As a first test, we demonstrate how the mixing and formation of 
secondary instabilities in the Kelvin Helmholtz instability (KHI) 
changes visibly from the original KHI tests that were presented in 
the Arepo paper (Springel 2010) using SSR. The setup of the ini¬ 
tial conditions of this shear fiow is described in Springel (2010). 
The evolved KHI at t = 2.0 is presented in Fig. 1 with both the 
LR and SSR schemes, at low and high resolutions. With the LR ap¬ 
proach, we have enforced near-Lagrangian motion with the correc¬ 
tion limiter parameter of / = 0.1. We ran the LR simulation with 
the old MH integrator and the GG gradient estimates. For com¬ 
parison, we ran the SSR simulation in the MH integrator and GG 
gradient estimates mode, as well as the new RK integrator and LSF 
gradient estimates mode. 

The SSR scheme shows evidence of mesh noise at both low 
and high resolutions; there is structure present on the scale of the 
mesh size, at either resolution. The LR method produces a result 
that is smooth and in agreement with high resolution fixed-grid 
simulations in terms of the secondary instabilities that develop (c.f. 
the static mesh runs shown for comparison in Fig. 1, or the static 
mesh runs in Fig. 34 of Springel (2010), or the static mesh finite- 
element simulations in Fig. 8 of Mocz et al. (2014)). Small scale 
secondary instabilities that are present with the SSR scheme disap¬ 
pear with the new LR method. The density field can be said to be 
well-resolved with the LR technique because there are no structures 
on the length scale of the cell size. It is worth pointing out that the 
solution is not exactly symmetric at high resolution because trun¬ 
cation errors lead to chaotic behaviour. We claim the secondary in- 
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Figure 2. Comparison of KHI with the SSR and LR schemes at low and high resolution. The LR scheme eliminates secondary numerical instabilities. 


stabilities present in the old SSR scheme runs are purely numerical 
artifacts, which we demonstrate below with a second KHI test. 

We investigate the KHI setup of McNally, Lyra & Passy 
(2012), which is a benchmark test in which only a single mode 
is excited and analytic theory exists to predict the growth rate. Very 
high resolution simulations of the instability are performed in Mc¬ 
Nally, Lyra & Passy (2012) using the 6-th order (in space) finite 
difference PENCIL code. In Fig. 2, we show that the SSR scheme 
produces secondary instabilities on the scale of the cell size, which 
are not present with the LR approach. The growth rate of the pri¬ 
mary mode of the instability of this test is still simulated accurately 


with both regularization methods and agrees with theory, as calcu¬ 
lated and presented in Fig. 3; only at late times (in the non-linear 
regime) are small deviations observed in the growth rate of the in¬ 
stability simulated with the two regularization techniques, where 
we expect the LR scheme is providing more accurate results due to 
the suppression of the secondary numerical instabilities. 

We explore the Lagrangian nature of the SSR and LR schemes 
in Fig. 4, where we plot the relative changes in the mass of each 
cell as a function of time. In a purely Lagrangian scheme, there 
is no mass exchange between particles. We see that the SSR and 
LR schemes saturate to the same level of mass exchange, regard- 
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Figure 3. Growth rates of the KHI simulated using the two regularization 
methods both show good agreement at early times with theory. 



Figure 4. Comparison of the mass of each cell as a function of time in a 
KHI simulation for the different regularization schemes (static mesh, SSR, 
LR with 1 /4 the time step, LR). The lines and shaded regions represent the 
quartiles of the distribution. 

less of whether the timestep is decreased by a factor of 4. Hence, 
we observe no strong implicit dependence of the deviation from 
Lagrangian behaviour on the time step. We have also explored re¬ 
ducing the number of Lloyd iterations from 5 to 2, and the mass 
exchange as function of time and the density solution at t = 2 do 
not change noticeably, so they are not shown. The mass exchange 
on a static mesh saturates to a higher level than the moving mesh 
approaches due to its Eulerian nature. 

3.2 Reduction of noise in the MHD Orszag-Tang Vortex 

We next consider simulations of the Orszag-Tang vortex (Orszag 
& Tang 1979), a common MHD test problem that involves super¬ 
sonic shocks and decaying turbulence. The setup of the problem 
is described in detail in Mocz, Vogelsberger & Hernquist (2014). 
We solve the system using the constrained transport algorithm 
for a moving mesh described in Mocz, Vogelsberger & Hernquist 
(2014). The outputs of the simulation at t = 0.5 are presented 
in Fig. 5. The SSR scheme produces numerical noise in the den¬ 


sity field, particularly at the locations of shocks and shear flows 
(compare the center and the shocks in the four cardinal directions 
in the image). The degree of these numerical noise errors can be 
changed by altering the mesh regularization parameters of the SSR 
technique, but it is difficult a priori to predict to what extent. For 
the tests presented here, we chose the typical values suggested in 
Springel (2010), and also explored stronger regularization parame¬ 
ters, but found that the noise does not go away completely because 
the corrections to the mesh generating point locations are not ap¬ 
plied in precisely the right directions to keep the cells centroidal to 
high precision. These numerical artifacts, however, are completely 
eliminated by the LR scheme, which produces a clean, smooth so¬ 
lution. We note that there are no parameters to adjust in the LR 
approach: the simulation was run without the optional Lagrangian- 
enforcer (i.e., / = oo). 

In Fig. 5 we tagged cells initially on the y = 0 axis at t = 0 
to show where they finish at the end of the simulation, to demon¬ 
strate that both regularization methods move the cells approxi¬ 
mately with the fluid flow, as they land approximately in the same 
places, demonstrating that the scheme is quasi-Lagrangian. Thus, 
just a tiny correction to the cell vertex velocities can have large 
implications on the amount of grid noise in the solution: directed 
at the right orientations they can be used to remove the grid noise 
almost entirely. 

We verify that the LR scheme keeps the mesh regularized to 
a much higher degree. A histogram of the relative center-of-mass, 
mesh generating point offset for the two regularization techniques 
is shown in Fig. 6 . The LR method shows over an order of mag¬ 
nitude improvement, and can be arbitrarily improved with more 
iterations in the regularization algorithm. The offsets in the SSR 
scheme, on the other hand, can be quite large (> 10 per cent of the 
cell’s effective radius). 


3.3 Proper 2nd order convergence in Yee Vortex and 
Improved Angular Momentum Conservation 

With the LR scheme, we can achieve formal second order conver¬ 
gence with a moving Voronoi mesh code, which has been long¬ 
standing issue (see also Pakmor et al. (2015) which achieves second 
order accuracy using an RK integrator and LSF gradient estimates). 
In many applications second-order convergence is not possible be¬ 
cause the astrophysical flows generate shocks, and in the presence 
of these shocks the slope limiter reduces the order of accuracy to 
first-order in order to maintain numerical stability across discon¬ 
tinuities. However, second-order convergence should be expected 
for smooth flows. 

We achieve second order convergence by using the LR scheme 
and the time-symmetric RK integrator (Pakmor et al. 2015) but un¬ 
like Pakmor et al. (2015) we can use GG gradient estimates since 
the mesh remains highly centroidal. The RK integrator averages a 
flux calculated with the mesh geometry at the beginning and end of 
the time step, instead of using just the mesh geometry at the begin¬ 
ning of the time step as the original MH integrator of Arepo does. 
We demonstrate the ability of the code with the new LR scheme to 
show second-order convergence in the Yee vortex problem, a 2D 
isentropic, differentially-rotating, steady-state smooth flow (Yee, 
Sandham & Djomehri 1999). 

The flow is described with parameters Tinf = 1, /3 = 5, 7 = 
1.4, and box size L — 10. The temperature, density, velocity, and 
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Figure 5. Density field at t = 0.5 in the Orszag-Tang vortex simulated with the SSR and LR schemes. The LR scheme eliminates noise at locations of shocks 
and shear flow in the simulation. A number of cells are tagged in black (initially on the y = 0 axis at f = 0) to visually verify that they have travelled to 
similar locations at the end of the simulation. 


pressure profiles are given as follows. 


T(r) = r„ - J 

(4) 

p{r) = 

(5) 

Vx{r) = 

(6) 

Vy{r) = 

(7) 

P[r) = h 

(8) 


We simulated this flow to t = 10, a time at which the vor¬ 
tex has differentially rotated several times. Fig. 7 shows a well be¬ 
haved profile using the LR scheme, but a noisy one with the SSR 
technique, even with the RK time integrator. We see that due to the 
presence of this mesh noise, the SSR method is unable to achieve 
second-order convergence of the LI norm error of the density pro¬ 
file, as shown in the convergence plot Fig. 8. Beyond some fairly 
low resolution, the mesh noise dominates the error in this test prob¬ 
lem. However, this noise is removed with the LR scheme. 

We simulate the flow much longer, up to t = 100 to study the 
long term behaviour of the angular momentum. The LI norm error 
of the angular momentum is reduced by an order of magnitude once 
the LR scheme is added, as shown in Fig. 9, and grows with time 
in a predictable manner as for this problem. The angular mo¬ 
mentum without the LR scheme can fluctuate unpredictably due to 
mesh noise, as also shown in Fig. 9. The total angular momentum 
is well behaved with the LR and RK time integrator improvements, 
as shown in Fig. 10, while with the original formulation of Arepo 


it can grow non-linearly, which has previously placed some limita¬ 
tions on the application of moving Voronoi mesh codes like Arepo 
or TESS for rotationally symmetric accretion flows. With the LR 
scheme, we have improved the angular momentum conservation of 
the moving mesh approach, and at the same time reduced the mesh 
noise. 

The well-behaved convergence with LR is achieved because 
there are no errors in the gradient estimate present that are caused 
by deviations from a fully centroidal mesh, which again result from 
the noise in the movement of the mesh-generating points. We em¬ 
phasise that the combination of GG gradients and a centroidal mesh 
even guarantees second order accurate gradients, whereas LSF gra¬ 
dients are only first order accurate in general (for a Cartesian mesh, 
they are also second order). 


3.4 Reduction of artificial power in turbulence test 

We present the results of 3D subsonic driven turbulence, whose 
setup is described in Bauer & Springel (2012), and was also pre¬ 
sented in Mocz et al. (2014), where the system was investigated 
using finite-element methods. The energy power spectrum for such 
a turbulent system is well described by a Kolmogorov power-law 
on the range of spatial scales that are resolved by the mesh. How¬ 
ever, the energy cascade is artificially limited by the sizes of the 
cells, and there is excess artificial accumulation of energy on these 
scales. This excess build-up typically has tended to be larger in 
moving mesh codes than in static mesh codes, with mesh noise be¬ 
ing a likely culprit as it adds fluctuations on the length scale of the 
cell size. In Fig. 11, we show that the artificial accumulation of en¬ 
ergy on the scale of the cells is reduced with the LR approach. The 
LR run used a Lagrangian enforcing parameter of / = 0.1, same 
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Figure 6. A histogram of the relative center-of-mass, mesh generating point 
offset for the two regularization methods. The LR scheme shows over an 
order of magnitude improvement. The errors in the SSR technique can be 
quite large, >10 per cent. 



Figure 8. Convergence analysis of the Yee vortex. Proper second-order con¬ 
vergence in the LI norm is achieved with the LR and the RK time integrator. 
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Figure 9. The angular momentum profile error as a function of time. The er¬ 
ror is reduced by an order of magnitude with the addition of the LR scheme, 
and grows in a predictable manner as 


Figure 7. The density profile of the Yee vortex at f = 10, with the SSR and 
time stepping scheme, the RK time stepping scheme, and the RK scheme 
with the LR method. The LR eliminates the mesh noise. 


as the KHI simulations, chosen so that the relative offsets between 
centres-of-mass and mesh generating points are reduced by at least 
an order of magnitude from the SSR scheme. 

We study the Lagrangian nature of the simulations by adding 
passive Monte Carlo tracers in the turbulence tests, as in Genel et al. 
(2013). The number of exchanges between cells, A^exch, that Monte 
Carlo tracers experience as a function of time, should be ideally 
zero in a fully-Lagrangian code. We plot the evolution of this quan¬ 
tity in Fig. 12. The LR and SSR methods show greatly reduced 
exchanges from a static mesh run, as expected. For a Lagrangian 
enforcing parameter of / = 0.1, the LR scheme exhibits somewhat 
more exchanges than the SSR scheme. The exchange number can 
be reduced by decreasing / to make the Lagrangian enforcement 
stronger, at the cost of weakening the strongly centroidal condition. 
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Figure 10. The total angular momentum is well-behaved with the RK time 
integrator and the LR scheme, but may grow exponentially with the SSR 
scheme. 
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Figure 11. Power-spectrum of driven subsonic turbulence for the two reg¬ 
ularization methods. On the smallest spatial scales (large k), the SSR tech¬ 
nique shows an excess build-up of power due to mesh noise, which is re¬ 
duced with the LR scheme. 



Figure 12. The number of exchanges between cells, A^exch? that Monte 
Carlo tracers experience as a function of time in the various turbulent box 
runs using different cell regularization schemes and at different resolutions. 
A lower number of exchanges indicates a simulation that is closer to fully- 
Lagrangian. Using a moving mesh instead of a static mesh reduces the num¬ 
ber of exchanges by over an order of magnitude. The LR method, in keeping 
the mesh centroidal to a high degree, is somewhat less Lagrangian than the 
SSR scheme. 


3.5 Adiabatic jet simulations with KHI 

As a test of the new regularization method in an astrophysical con¬ 
text where it is expected to have an impact, we present simulations 
of heavy, supersonic, adiabatic jets that have a helical mode ex¬ 
cited which triggers the KHI, following the setup of Bodo et al. 
(1998); Micono et al. (2000). This situation is ideal for the study of 
different phases of the temporal evolution as a function of sound¬ 
crossing time. The jet has a surrounding gas density to jet density 
ratio of V — 0.1, and a supersonic Mach number of M = 10. 
Our simulations are in 3D and use 5 million cells. In Figure 13, we 
present the results of the simulation at a point where the KHI has 
developed and shredded the jet. The LR scheme shows lower level 
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Figure 13. Simulation of an M = 10, heavy, adiabatic jet with the two 
regularization methods. The cross-sections are at time t = 7 sound-crossing 
times. Entropy generation is higher at late times with the SSR scheme than 
with the LR scheme. 


of mesh noise and entropy generation than the SSR technique at 
late times in the simulation. 


3.6 Isolated galaxy formation 

Finally, we simulate the formation of an isolated disk galaxy with 
magnetic fields similar to the Milky Way, following the setup of 
Pakmor & Springel (2013). The initial gas sphere in the simulation 
is in hydrostatic equilibrium (without radiation), but collapses due 
to radiative cooling. The gas has some initial angular momentum, 
and settles into a dense, rotationally supported disk. Regions in the 
disc can fragment and are allowed to form stars. In this simulation 
with collapse, we bias the initial velocities of the mesh generating 
points by the density gradient term as in Vogelsberger et al. (2012), 
before applying the Lloyd regularization, in order to closely follow 
the collapse. We present the results of the simulations under the LR 
and SSR schemes as well as time evolution of some global prop¬ 
erties in Figure 14, and find the results are unaffected by the regu¬ 
larization because mesh noise errors are sub-dominant. The simu¬ 
lations have a mass resolution of 10® Mq. 
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Figure 14. Simulation of the formation of an isolated magnetic disc using the two regularization schemes. The first two panels show a slice of the gas density 
held, while the third panel shows the time evolution of mass of stars and gas within a 15 kpc radius. Regularization does not affect the global properties of the 
simulation. The density scale covers — 5 < log^^Q P ^ —0.8. 


4 CONCLUDING REMARKS 

In summary, we have presented a new regularization scheme, Lloyd 
Regularization (LR) for moving Voronoi mesh codes that signifi¬ 
cantly reduces mesh noise and improves convergence and angular 
momentum conservation. These advantages are gained because the 
new LR scheme keeps the cells centroidal to a high degree by ap¬ 
plying a forward-predicting Lloyd iteration correction to the mesh 
generating point. The regularization scheme will improve the sci¬ 
ence applications of moving mesh codes used in astrophysics, par¬ 
ticularly in systems where shear flows exist. The LR scheme allows 
moving mesh codes to capture the physics of fluid instabilities more 
accurately than earlier methods, which are thought to be important 
for jet dynamics and gas stripping in galaxies. The correct level 
of turbulence and entropy generation is also improved with the re¬ 
duction of mesh noise through this method. The improvements are 
most relevant for shear flows, and complements the recent moving- 
mesh improvements for smooth solutions of Pakmor et al. (2015). 
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